Lissajous curves as aerial search patterns

Manned and unmanned systems are prevalent in a wide range of aerial searching applications. For aircraft whose trajectory is not or cannot be planned on-the-fly, optimal deterministic search pattern generation is a critical area of research. Lissajous curves have recently caught attention as excellent candidates for all kinds of aerial search applications, but little fundamental research has been done to understand how best to design Lissajous pattern (LP)s for this use. This paper examines the optimization of these search patterns from analytical, numerical, and data-driven perspectives to establish the state of the field in Lissajous curves for aerial search. From an analytical perspective, it was found that the average expected distance between a Lissajous searcher and a random target on a unit square approaches 0.586 as search time increases. Furthermore, an analytical approximation for the average searcher speed was found to guarantee error of no more than 22.1%. Important outcomes from the numerical optimization of Lissajous search patterns include the development of an intuitive evaluation criterion and the conclusion that irrational frequency ratios near 0.8 typically yield highest performance. Finally, while a robust predictive model for fast pattern optimization is yet out of reach, initial results indicate that such an approach shows promise.


Background and motivation
The use of unmanned aerial systems (UAS) continues to proliferate, and thousands of new applications surface each year.The low cost of operating most UAS and the low barrier to entry has enabled operators to find uses for autonomous flying vehicles that are still relatively new and potentially disruptive to both commercial and government sectors.According to the International Civil Aviation Organization, there were over 2 million extant UAS in 2021 and this number is expected to rise to over 6.5 million in 2030 1 .UAS have been successfully incorporated into a host of activities including agriculture, aerial photography and videography, construction, mining, mapping and surveying, delivery services, and emergency response 2 .

Related work
Using UAS to assist in search and rescue has been an active area of research for over a decade.A key challenge for any autonomous or semi-autonomous system, especially in the context of aerial search, is path planning.Path planning can occur before or during missions, and can be updated in response to new information gained during the mission.Such information theoretic approaches to path planning are an area of active research 12,13 , and although adapting a path in light of information gained during search sounds attractive, there can be several practical drawbacks.One of the most significant hindrances in aerial search contexts is the excessive computational cost of continuously updating the search pattern, especially for small UAS with limited onboard computational assets and power.Moreover these patterns can be ineffective when little information about the search domain is known or made available.
Many recent ideas have been proposed in the areas of path planning as well as maintaining control of multiple unmanned aerial vehicles (UAV) flying simultaneously in support of a search mission.In 2009, Lin and Goodrich proposed an intelligent path planning method for wilderness search and rescue, however their technique relies heavily on a pre-determined probability distribution map to drive the performance of their search 14 .More recently Baker et al. 15 introduced a coordinated Monte Carlo tree search algorithm which showed good performance against more basic search methods, but the method relies on knowledge of the situation on the ground prior to the start of the search.Although appropriate for post-disaster scenarios such as after an earthquake event, it is less well-suited to more dynamic scenarios.
Hayat et al. proposed a multi-objective optimization algorithm to allocate tasks and plan paths for a team of UAVs 16 .This work envisions a large search area where communications can only be assured through a mesh network of air vehicles which also utilize sensors to aid in their search.The flight paths are updated to maintain network integrity, but the path itself is not updated based on information from onboard sensors.San Juan and colleagues 17 introduce a discrete path planning technique using four separate strategies, all of which assumes some level of prior knowledge about the location of the subject of the search based on a Risk/Occupancy Map, terrain, or other data.Research by Rahmes et al. 13 focuses on creating a probability map which is updated at each last known position to calculate the next position for the UAS which maximizes the likelihood of locating the target in order to avoid UAS spending too much time searching in low-probability areas.
Echeveste et al. 12 propose a UAS planning method with a focus on mapping the concentration of ground contamination by using Kriging variance to estimate the concentration in the entire area of interest based on five random sample points.In order to achieve the greatest reduction in uncertainty with each additional sample point, the study uses Variance Driven Sampling (VDS) to sample subsequent points with the greatest variance.However, this method is computationally expensive and may require the UAS to travel the across the entire search space to reach its next sample point.
Finally Xing et al. 18 developed a multi-UAV cooperative system for search and rescue based on the 5th Generation You Only Look Once (YOLO) algorithm.Their framework includes the ability to free-graft multiple UAVs, independent control of each UAV, real-time target detection, and monocular positioning.Although sophisticated and able to accommodate heterogeneous search agents, the algorithm relies on hand selected initial flight paths and does not perform real-time updates to the path itself in spite of agent-to-agent communication with search status updates.There remains a lack of investigation into deterministic patterns that are both effective and efficient at finding targets in different scenarios.
Historically employed deterministic patterns are outlined by the International Aeronautical and Maritime Search and Rescue Manual 3 .Figure 1 shows examples of the most common patterns that have proven to be effective for maritime applications.In general, deterministic search patterns can often be classified as either cyclical or comprehensive.Cyclical patterns are usually designed for contexts where a moving target might be intercepted, and as such, they commonly do not cover an entire search space.In contrast, comprehensive patterns are designed to methodically search a whole area to guarantee the target will be located by the completion of the search.A weakness of comprehensive patterns is that they trade a potentially faster detection time for higher confidence of overall mission success.The hybrid pattern proposed in 19 uses Lissajous curves to execute aerial search.This approach combines cross-domain sweeping motions with methodical area coverage to provide a thorough yet efficient interrogation of a search area.Search methods based on Lissajous curves can be both cyclical and comprehensive as demonstrated in Fig. 2, and they show significant promise as a potential alternative search method to those in current use 20 .In particular, LPs have been shown to more quickly reach an 80% confidence threshold of mission completion when optimized for assumed target size 21 and are especially effective against canonical patterns when the target is modeled as adversarial 22 .Lissajous flight patterns can be quickly generated, are easily tuned to accommodate a wide range of scenarios, and do not require oversight or input about target location probabilities.

Original contributions
The primary objective of this research is to comprehensively investigate the design of Lissajous search patterns for maximum efficacy from analytical, numerical, and data-driven perspectives.Original contributions include not only insights into optimal LP design for UAS search, but also a few mathematical derivations that are more generally useful in their own right.These include: (1) generalized continuous and discrete formulations for a Lissajous curve arbitrarily located and rotated in two dimensions, (2) the expected distance between a known searcher location and a randomly located target uniformly distributed on a unit square, (3) a symmetric linear approximation of the Pythagorean theorem, and (4) an approximate relationship between the average speed of a Lissajous curve and the parameters that describe it.Section "Lissajous curves" first establishes what Lissajous curves are, including a few existing formulations and some discussion about their features.Everything that follows is put forward as novel work: section "Analysis of LPS for aerial search" addresses LPs from a purely analytical perspective, motivating the need for numerical methods which are detailed in section "Numerical assessment and design of Lissajous search patterns".Section "Predictive modeling for Lissajous search pattern design" investigates the use of predictive models for more efficient hyperparametric pattern design, and section "Conclusions and future applications" conveys concluding thoughts and future directions of the research.

Traditional search patterns
While aerial vehicles can be controlled in a variety of ways, waypoint navigation is perhaps the most flexible and widely used.Waypoint navigation requires the generation of discrete locations in space which an agent must visit at specific times.The density of waypoints required to sufficiently construct a search pattern depends on the complexity and curvature of the pattern.The following subsections describe waypoint generation for the four canonical search patterns used almost exclusively in the field today, depicted in Fig. 1.It should be noted that all the patterns discussed here assume a constant searcher altitude.While three-dimensional search certainly has importance, we leave 3D considerations for future work.

Circle search
Because circle search (CS) patterns can only be well approximated by many straight line segments, they may require well over ten waypoints to adequately construct.The x and y coordinates for a circular search agent are obtained by the following: where R is the radius of the search circle, v is the velocity of the agent, [x ȳ] ⊤ are the coordinates of the center of the circle, and t is the time the searcher takes to move between waypoints.This time step parameter controls the resolution of the pattern.

Sector search
The sector search (SS) path uses a cyclical, overlapping pattern whose waypoints are most easily defined by a recurrence relation.The k th x (east) and y (north) coordinates and two-dimensonal (2D) heading θ of a search agent are given by: where and D is a distance determined by the breadth of the search region.This pattern repeats when k is a multiple of 6, at which point an angular offset is often added.

Parallel track
The parallel track search (PTS) pattern methodically covers a rectangular search space by scanning the width or height of the space at regular offsets.The waypoints for a PTS moving from south to north can be generated as follows: where (1)

Expanding square
The expanding square (ES) begins at or near the center of a search domain and moves outward along a spiraling path.The critical waypoints for such a pattern are given by: where and l is a distance parameter governed by the size of the search space.

Continuous parametric formulation
Despite their potential for complexity, Lissajous curves have a relatively simple mathematical formulation in two dimensions.Equation (8) give the parametric representation of these curves in the x and y directions: where A x,y are the amplitudes (half-widths) of the pattern in the horizontal and vertical directions, ω x,y are the angular frequencies of the pattern, and φ x,y are the phase shifts which dictate the starting coordinates of the curve.By tuning these six parameters, any rectangular Lissajous curve can be generated, though this formulation restricts the pattern to be centered at the origin and orthonormal with the x and y basis vectors.

Discrete recursive formulation
A simple version of the discrete linear recurrence relation describing 2D Lissajous curve waypoint generation was derived in 19 .This linear discrete recursive formula is a useful representation of the curve which is amenable to linear estimation and control frameworks.Define the discrete state vector x k as follows: This state is propagated according to the following linear model with sample period t: where The matrix A is a second-order approximation of the state transition matrix which precisely discretizes the Lissajous curve.The desired amplitudes and phases of the pattern are used to derive the initial state conditions:

Features of Lissajous curves
In addition to the formulations above, a qualitative description of Lissajous curves is beneficial to better convey their utility for aerial search.Figures 3-6 are included here to visually demonstrate the effect of varying each of the key parameters in Eq. ( 8).As Fig. 3 demonstrates, the amplitude parameters A x and A y of Eq. ( 8) can be adjusted to match the width and length of a rectangular search space.
While the effect of the amplitude parameters on the shape of the curve is relatively straight-forward, the effect of the phase shifts φ x and φ y is more complex.These parameters determine the start point of the pattern and they also have an influence on the size of gaps in the curve.This is demonstrated by Fig. 4.
When the area of the smallest gap collapses to zero, the pattern takes on a lower order of complexity (this can be seen in Fig. 3, where φ x approaches zero and φ y approaches π/2 ).Furthermore, as Fig. 4 shows, a given (6) combination of φ x and φ y does not uniquely specify a Lissajous curve.For a fixed set of amplitudes and angular frequencies, there are infinitely many φ x,y combinations that yield the same pattern, but there are infinitely more phase shifts which produce different curves.However, the number of intersections in a Lissajous curve (i.e. its level of complexity) can only be controlled by varying the frequency parameters ω x and ω y .More specifically, it is the ratio of the angular frequencies in the x and y directions that governs complexity, as the absolute value of these parameters affects only the "speed" of the pattern.Figure 5 demonstrates examples of three different frequency ratios.
It is worth noting that frequency ratios near one cover a space with higher rotational symmetry than smaller ratios: this is evident when comparing the 2/7 pattern with the 24/25 pattern.Finally, it can be shown without   www.nature.com/scientificreports/much difficulty that irrational frequency ratios result in curves which never repeat.Figure 6 contrasts an irrational frequency ratio with a nearly equal rational one.

Formula generalization
The first original contribution of this work is a generalization of the LP which allows the curve to be rotated and shifted in 2D space.

Continuous
The generalized continuous formula for a Lissajous curve is given below: where X and Y are the offsets of the pattern in the x and y directions, respectively, and θ is the angle by which the bounding rectangle of the pattern is oriented.This formulation offers more versatility in constructing Lissajous curves than Eq. ( 8).

Discrete
Rotating and translating a Lissajous curve using the discrete recursive representation is more challenging than in the continuous case.While a more thorough derivation can be found in "Appendix A", the result is given here.Defining the state vector as in Eq. ( 9), the state is propagated by the following discrete linear model: (11), and This also requires that the initial state be transformed according to: This construction of Lissajous curves allows any pattern to be generated at any location and orientation in the plane.

Expected distance to a uniform random target
Designing an "optimal" Lissajous search pattern is an open-ended task.One must take into account many factors, including the assumed prior belief about the target, the dimensions of the search space, the maximum duration of the search, and so on.To begin to approach this problem analytically, several assumptions must be made about the circumstances of the search problem at hand.The developments of this section make the following assumptions about the search scenario: 1.No prior knowledge is available about the target's location, except that... 2. ...its distribution has 90 • rotational symmetry, and...For generality, all variables in the following developments are dimensionless.To further simplify, let θ = 0 and X = Y = 0 without loss of generality.Assumption 3 is inherently compatible with a Lissajous curve, as these patterns are always inscribed by a rectangle.From assumption 4, the values of A x and A y in Eq. ( 13) must be 0.5.In order to enforce assumption 5, x 0 must be [0.5 0 − 0.5 0] ⊤ .From this, it follows that φ x ≡ π 2 (4n − 3) and φ y ≡ − π 2 (4n − 3) , or for simplicity, φ x ≡ π 2 and φ y ≡ − π 2 .Finally, to simplify the problem formulation in keeping with assumption 2, let the frequency ratio r ω ∈ R * + be introduced and defined as follows: This is the single most critical parameter in varying the complexity of an LP, as described in section "Features of Lissajous curves".The assumption that the target distribution has 90 • rotational symmetry allows the range of values taken by r ω to be truncated to the half-open interval (0, 1].Thus, reflecting assumptions 2-5, the continuous Lissajous formulation of Eq. ( 13) can be reduced to the following: (Note: the S subscripts on X and Y denote the x and y positions of the searcher.)Therefore, with only a few reasonable and nearly trivial assumptions, a nine-dimensional parametric optimization space is reduced to only two.While assumptions 2-5 affect the parameters of the LP used in a search scenario, the first assumption relates to the assumed target distribution.If no information is available apart from the boundaries of the search domain, the distribution chosen to probabilistically represent the location of a target must maximize the differential entropy of the target.It can be shown 23 , and is also fairly intuitive, that the uniform distribution accomplishes this.If the support of the target location is known, the boundaries of the search pattern should be made to match in order to maximize the efficiency of the search.Since assumption 4 restricts the dimensions of the Lissajous curve, the probability density function (PDF) of the target's location X T , Y T must be: Analytically optimizing an LP subject to the above assumptions and constraints becomes a matter of determining which value(s) of r ω result in the fastest completion of a search mission.Without inserting any additional information about the modality of target detection, it may be reasonably assumed that such an optimal pattern should be made to minimize the expected distance between the searcher and the target, averaged over all time.The distance between the searcher and the target in two dimensions is simply (altitude is removed from the problem since it is assumed that the searcher can only identify the target when it is nearly directly above it).The expected value of this distance is: This is a notoriously challenging integral.However, it may be noted that the integrand is more easily represented in polar coordinates, in which case the integral becomes: where r = (x S − x T ) 2 + (y S − y T ) 2 .While the operation is now simpler, the region of integration is no longer constant.It is useful to visualize the surface being integrated to gain insight about the integration limits; this visualization is provided in Fig. 7.
It can be shown that the integral of the white-highlighted triangular region of the figure is: . This expected distance function is minimized to a value of 0.3826 when x = y = 0 , maximized to a value of 0.7639 when x = y = 1 2 , and is rotationally symmetric.The surface described by Eq. ( 24) is plotted in Fig. 8.This figure also depicts the surface visually decomposed into the α and β terms alone (It may be helpful to note that α and β are not significant functions themselves but are just defined here for compactness.They are plotted in the figure to show the relative contribution of each to the total expected distance.).
Let D(r ω ) represent the average expected distance from a searcher to a random target over a long period of time:  If a value of r ω can be found which minimizes D(r ω ) , it may follow that an efficient Lissajous search pattern should be designed with that minimizing frequency ratio, as it guarantees that the searcher will be as near to the target (on average) as possible.However, combining Eqs. ( 13) and ( 24) is extremely unwieldy and the integral of Eq. ( 25) is intractable.Nevertheless, it is possible to evaluate D(r ω ) numerically.The average expected distance was computed over 10,000 irrational values of r ω from 0 to 1 to ensure that the pattern never repeated.This was conducted for five values of T ranging from 10( 2π ω y ) to 100( 2π ω y ) as shown in Fig. 9 and summarized in Table 1.Though it cannot be proven analytically, it is clear from the figure that D(r ω ) converges to a value of approximately 0.586 as T → ∞ .This is further supported by Table 1.Therefore, we conjecture that over a long time, the average distance between a Lissajous searcher and any target uniformly randomly located within a square region is the same regardless of the shape of the LP.This statement is further supported by the fact that the average expected squared distance can be proven to be a constant 5 12 on the support of the unit square (see "Appendix C").While it is an interesting and somewhat unintuitive result that the average expected distance to a randomly located target is independent of the searcher's Lissajous path, it does not aid in the design of an optimal Lissajous search pattern.Furthermore, this approach cannot account for the technicalities of target detection (i.e. at what distance between a target and a searcher can the target be considered found?).For these reasons and others, practical Lissajous search pattern design must be done numerically.Before this is detailed in section "Numerical assessment and design of Lissajous search patterns", it is critical to know the relationship between an LP's parameters and its average speed in order to conduct fair and consistent numerical simulations.This can be thought of as a normalization of Lissajous patterns.For example, if numerical optimization shows that low values of r ω appear to perform the best, yet those frequency ratios cause a searcher to travel much faster, on average, than other values of r ω , the comparison would be misleading.Therefore, a derivation of the average speed of a Lissajous searcher is first appropriate.

Average speed of a Lissajous searcher
The instantaneous speed (i.e. the magnitude of the velocity vector) of a continuously formulated LP subject to the assumptions of section "Expected distance to a uniform random target" is given as follows: where dots denote time derivatives.The average speed over some period T is then   www.nature.com/scientificreports/This integral is again intractable.However, approximations can be considered to varying degrees of accuracy.The alpha-max-plus-beta-min algorithm commonly used to speed up magnitude computations in digital signal processing [24][25][26] offers the following approximation for s(t): where optimal values of α = 0.96 and β = 0.40 guarantee an error of less than 4%.However, this algorithm requires that the larger of ẋ(t) and ẏ(t) be known for all t-knowledge which cannot be guaranteed for a general Lissajous curve.Thus, a symmetric approximation is sought."Appendix D" shows that the optimal symmetric linear approximation of the Pythagoream theorem is: with error no greater than 22.1% .While it introduces a larger error margin, this approximation enables a much simpler estimate of the average speed of a Lissajous searcher.
The estimated average speed s is determined by Let the interval T over which the speed is being averaged allow for a complete number of periods in both the x and y directions: where m ∈ Z is the number of cycles in the x and n ∈ Z is the number of cycles in the y (While these developments are technically only valid for rational frequency ratios r ω = m n , the errors induced when applying the result to irrational frequency ratios are negligible.).It can be shown without much difficulty that and so, Solving for ω y (with ω x following from Eq. ( 17)) gives (More generally, for an LP with arbitrary amplitudes, Therefore, for a desired average Lissajous searcher speed of s and a frequency ratio of r ω , it can be guaranteed that choosing ω y and ω x using Eqs.(34) will result in a true average searcher speed with an error of no more than ≈ 22% .To confirm this, 40,000 LPs were generated with random values of ω y ranging from 1 to 100 and frequency ratios ranging from 0 to 1. Equation (33) was applied to solve for the estimated path speed and this was compared to the "true" (numerical) measure.The percent-error surface is plotted in Fig. 10.Note that the upper limit of the error indeed approaches the theoretical limit as ω y → 0.

Monte Carlo simulation
Given the analytical complexity of designing LPs, a Monte Carlo simulation technique was implemented to numerically achieve Lissajous flight path optimization.As described above, the goal of LP design is to choose a frequency ratio r ω that produces the most effective search path subject to the assumptions and constraints of section "Expected distance to a uniform random target".However, with a numerical/simulation approach, one can bypass the task's analytical intractability to deliver waypoints for aerial vehicles conducting real-world search missions.Multiple aerial searching scenarios can be simulated simultaneously by populating a unit square space with uniformly randomly generated targets, represented by their 2D Cartesian coordinates.Previous work 22 showed that 10,000 random targets are generally sufficient to achieve convergence of outcomes for the constraints of the scenarios addressed here.The Lissajous search path is then numerically represented by a sequence of dense waypoints generated by Eqs.(13) or (14).To ensure equitable search pattern comparison, waypoints were generated at resolutions that delivered constant average Lissajous speeds by calculating the appropriate values of ω x and ω y for each pattern.Although an approximate analytical approach for this was derived in section "Average  (31) speed of a Lissajous searcher", a pre-processed numerical lookup table was instead used here to reduce potential maximum error from near 22% to only fractions of a percent.At each waypoint, the detection state of each target was evaluated.Unlike in the analytical case, it is possible in simulation to account for the modality of target detection when assessing the performance of a search pattern.The simulated detector checks which targets lie in a rectangular field of view (FOV) at each waypoint in order to emulate the detection of a small target by an optical sensor.To limit the number of additional parameters introduced in the optimization problem, we assume a square FOV as projected onto the ground plane with side length l f .Detection is then accomplished efficiently in a two-step process.First, the Euclidean distance between the searcher and all targets is computed and targets whose distance from the searcher exceeds √ 2l f (the radius of the circle which circumscribes the FOV) are removed from further consideration.The much smaller subset of targets which lie near the agent at a given waypoint are determined to be in the FOV by first transforming their coordinates into the searcher's reference frame.If the magnitude of any of these transformed coordinates is less than l f /2 , the corresponding is considered to have been detected.
Figure 11 illustrates such a search scenario with uniform random target locations, searcher path, and FOV displayed.As one might expect, there is an important interplay between the maximum gap size in an LP and the size of the FOV.This relationship, investigated more thoroughly in the results of 22 , is taken into account in the statistical modeling of section "Predictive modeling for Lissajous search pattern design".
Even when a target lies in a UAS's FOV, no real-world detector can deliver 100% certainty of finding it, as sensor-based detection is always ultimately stochastic.For this reason, a probabilistic model for detection was chosen based on 27 to transform the time t d a target spends in the detectable region (i.e.FOV) of a searcher into  www.nature.com/scientificreports/ a probability of having been detected.The longer a particular target is inside the FOV, the more likely it is to be seen by the simulated UAS.This probability is modeled as where τ is a scaling time constant which holistically models the sensitivity of the detector (a smaller time con- stant corresponds to a more robust detector).At the conclusion of a simulation, each target has an associated probability of having been detected by the end of a search.Additional insight can be obtained by generating the cumulative density function (CDF) F T (t) which describes the probability that the time at which a uniformly randomly located target is found (detected) at a time T is less than t.Examples of Lissajous search pattern CDFs are shown in Fig. 12.While this figure is intended to serve primarily as a qualitative demonstration of how the performance of search patterns affects the shape of their CDFs, an accompanying quantification of this figure is provided in Table 2, which gives the CDFs' means and final values.
The CDF is a helpful visual representation of a particular LP's performance in the time domain, but it does not provide a single summative value to use as an objective function in the optimization problem.The mean or final value of the CDF may be reasonable evaluation metric candidates, but neither of these parameters captures much information about the early behavior of the pattern.For this reason, an evaluation criterion which reflects the insightful shape of the CDF is desired for comparison between LPs.

Evaluation criteria
Previous work 21 proposed the use of the area under the CDF (AUC) as a qualitative summary metric for characterizing the effectiveness of a search pattern.If a pattern is particularly robust, its CDF will demonstrate a steep rise early on, contributing to a large overall AUC.However, the challenge introduced by this metric is the arbitrariness of the upper limit of integration.It was formerly proposed that this limit should be chosen as the time  at which a certain acceptable threshold for the probability of detecting a target has been achieved.Maximizing the AUC has practical utility in providing an optimization objective, but it lacks a clear underlying mathematical interpretation relating to the stochastic problem at hand.In this work, we instead propose to use the area above the CDF (AAC) as an objective function.Maximizing the AUC and minimizing the AAC achieve the same goal, but using the AAC has two distinct advantages over using the AUC: (1) the arbitrariness of a cutoff for integration is eliminated since a CDF always converges to one (This is always true if there are no gaps when the FOV is swept out over the search pattern.), and (2) the area above a CDF which is defined strictly over a non-negative domain is in fact equal to the mean of the distribution.This gives meaning and mathematical significance to the AAC as a criterion for evaluating a Lissajous search pattern.
The area above the CDF, or the mean of the distribution, can be shown 28 to be: where E[T] is the time at which a uniformly randomly located target can be expected to be found by a UAS searcher flying on a prescribed Lissajous path.Therefore, smaller E[T] values are associated with more effective frequency ratios for pattern design.It is noteworthy that the upper integral bound in Eq. (36) need not necessarily be ∞ .If a mission highly prioritizes the speed with which a target is located over the certainty of finding the target at all, a lower certainty threshold may be acceptable in which case the upper bound of the integral would change.If the desired certainty level for mission success is α , Eq. (36) becomes where now E[T α ] can be described as "the expected time at which a search mission for a uniformly randomly located target is successful".Figure 13 demonstrates this concept for α = 0.85 , a typical certainty threshold for urgent search missions.Previous work has shown that Lissajous curves are particularly advantageous over traditional deterministic patterns in such scenarios when 100% certainty of locating the target is not required 22 .

Key findings
Monte Carlo simulations were conducted as described above.Given the assumed 90 • rotational symmetry of the assumed target distribution, only frequency ratios spanning the interval (0, 1] require investigation, thus increasing the resolution of numerical optimization without introducing unnecessary computational burden.Figure 14 shows the dependency of LP performance on frequency ratio for α = 0.85.The absolute minimum E[T α ] over this range naturally corresponds to the most effective frequency ratio for an LP, where any local minimum is a well-designed pattern and any local maximum is poorly suited for aerial search.One key finding from this approach is the recurring poor performance of Lissajous paths with near-rational frequency ratios.The large spikes at rational intervals imply that a higher degree of irrationality in frequency ratio is indicative of higher-performing LPs.This is an intuitive result in light of the fact that irrational patterns will never repeat, so "more rational" patterns will wastefully overlap themselves more often.Furthermore, rational frequencies are prone to leaving more unsearched areas in the search space as shown in Fig. 6, thus resulting in higher mission completion times.With these results, the task of designing a Lissajous search path becomes simply selecting the frequency ratio corresponding to the optimal minimum E[T α ] value.This pattern will be the most effective flight path given the parameters of the search scenario that was simulated.
While the simulation method yields important design results, it is not without drawbacks.The numerical optimization process is computationally intensive, discouraging real-world implementation where deployment . A CDF evaluation where the mission's desired certainty threshold is 85% rather than 100%.
speed is of high priority.For this reason, a predictive modeling technique was considered for frequency ratio selection which bears the bulk of the computational cost up-front at model fitting, but almost entirely eliminates the burden at the implementation stage.

Data and methodology
As Fig. 14 shows, under most conditions there are several ranges of nearly-optimal frequency ratios for which search performance varies only minimally.In fact, most values of r ω result in moderately acceptable performance as compared with the assorted outlier "spikes".Thus, an alternate approach to Lissajous search pattern design may ask not the question, "What frequency ratio is optimal?" but rather "What range(s) of frequency ratios should be avoided?"When a rigorously exact globally optimal solution is not required, data-driven approaches can provide an adequate solution much faster than pure simulation alone.This is afforded by the fact that such a predictive modeling approach bears the computational burden of simulation in offline training during model generation rather than online querying.The frequency ratio r ω is an explanatory feature which only partially explains the response variable E[T α ] , but there are several other features (subject to the assumptions of section "Expected distance to a uniform random target") that also affect the response.In particular, these include: (1) the size of the agent's FOV as a percentage of the area of the search space, (2) the time constant τ from Eq. 35, and (3) the certainty threshold α from Eq. 37.
To build the predictive models, 86,999 simulations were executed with each explanatory variable being randomly generated and the response variable recorded.Ten percent of these trials were retained as a testing set, while the rest were used to train predictive models.Table 3 summarizes the statistics of the simulations.The range of values for %FOV, α , and τ were chosen to roughly represent realistic scenarios.

Modeling
To begin, a basic linear regression was used to predict the expected mission completion time from the four simulation input variables.After building and testing the linear model, only 35.74% ( R 2 = 0.3574 ) of the variation in the data was explained, and the root-mean-squared-error (RMSE) in predicting the test set was 1573.542 as shown in Table 4. Including the interaction terms between each of the four explanatory features only marginally improved model performance, increasing the R 2 to 0.3929 and decreasing the RMSE to 1519.414, still a poor model.However, linear regression inherently assumes a normal distribution of the response, and as Fig. 15 shows, the distribution of E[T α ] is much more resemblant of a gamma distribution.
Taking this clue, the linear model was generalized to use a logarithmic link function and assume an underlying gamma distribution of the response.The resulting gamma regression model provided improved results with  a pseudo-R 2 ( R 2 L ) of 0.5609 and RMSE of 1437.222.The R 2 L measure relies on a likelihood ratio, hence the "L" subscript, and is calculated by where D 0 is the model's null deviance and D k is the residual deviance with k degrees of freedom 29 .
While gamma modeling is the most appropriate regression generalization for the problem at hand, basic regression is still fundamentally insufficient to adequately model the Lissajous simulation data.Traditional gamma modeling with a logarithmic link function assumes a linear relationship between the explanatory features and the logarithm of the response, but Fig. 16 shows clear evidence of a nonlinear relationship.
However, a generalized additive model (GAM) provides a more flexible method that can characterize nonlinear regression effects 30 .While a basic linear model assigns to each predictor an associated coefficient β i : a generalized additive model fits a non-parametric function f i (natural cubic splines in this case) to each predictor as shown in Eq. (40): A GAM allows the modeling of nonlinear effects while still leveraging an underlying gamma distribution to achieve better results than either previous technique.A GAM using natural splines has an additional hyperparameter for each non-parametric function it fits, realized as the number of "knots" for each spline.Using five knots for each function yields an R 2 L of 0.7052 and an RMSE of 1270.333.Increasing the number of knots can increase R 2 L , apparently improving the model, but this also increases the RMSE, which is a sign of over-fitting.Thus, five knots were heuristically chosen for the results delivered here.A summary of the applied statistial models is given in Table 4.

Figure 1 .
Figure 1.Four common deterministic search patterns prescribed in the International Aeronautical and Maritime Search and Rescue Manual 3 .From upper-left to lower-right: circle search, sector search, parallel track search, and expanding square search.

Figure 2 .
Figure 2. Assorted Lissajous curves.See section "Features of Lissajous curves" for a more detailed explanation.
14:11144 | https://doi.org/10.1038/s41598-024-60803-2www.nature.com/scientificreports/This is worked out in "Appendix B".Adding the result of this integration with those of the remaining three regions gives the following result for the expected distance E[d(x S , y S )] between a uniformly randomly distributed target location on the unit square and a given searcher location (x S , y S ):

( 24 )Figure 7 .
Figure 7. Surface of integration for distance function d(x S , y S , x T , y T ).

Figure 8 .
Figure 8. Expected distance between searcher and randomly located target.

Figure 11 .
Figure 11.Simulated UAS with square FOV navigating the search space along a r ω = 0.3 Lissajous path.

Figure 12 .
Figure 12.Monte Carlo simulation results showing examples of search CDFs and their corresponding Lissajous curves over a range frequency ratios.Curves should be read in increasing order from upper-left to lower-right.For these simulations, τ = 20.

Figure 14 .
Figure 14.Frequency ratio decision plot showing the expected time to complete a search mission for every simulated LP (1000 frequency ratios tested).Optimal r ω shown at ∼ 0.779 with E[T α ] = 564.926.

Figure 15 .
Figure 15.Distribution of the independent response variable E[T α ].

Figure 16 .
Figure 16.Strongly nonlinear relationship between frequency ratio r ω and the logarithm of the response E[T α ].
odd , W is the width of the search region, and H is the desired distance increment in the y direction.

Table 2 .
Mean and final value of CDFs for example Lissajous patterns.

Table 3 .
Training set summary statistics.